Field survey data

Using the package glmmTMB to model groups of species individually.

The first 6 rows of the dataset
month site salinity_ppt salinity_regime quadrat_no balanus_pt chthamalus_pt barnacle_recruits_pt barnacles_total_pt mytilus_pt fucus_pt ulva_pt hildenbrandia_pt mastocarpus_pt mastocarpus_crust_pt pyropia_pt endocladia_pt microcladia_pt polysiphonia_pt diatoms_pt urospora_pt acrosiphonia_pt melanosiphon_pt leathesia_pt greens_pt greens_pyropia_pt reds_pt reds_fucus_pt littorina_spp_no sitkana_no littorina_total_no paradigitalis_no unknown_limpet_no pelta_no scutum_no persona_no digitalis_no all_limpets_no gastropods_no anemone_pt hemigrapsis_no hermit_crab_no
May Ruckle 30 H 1 1.00 0.50 0.0 1.50 0.250 0 0.00 0.00 0.00 0.50 0.00 0 0.00 0 0 0 0 0 0 0.00 0.00 0.50 0.50 41 0 41 49 0 0 0 0 2 51 92 0 0 0
May Ruckle 30 H 2 8.25 13.00 0.0 21.25 0.025 0 0.00 0.50 0.00 0.00 0.00 0 0.00 0 0 0 0 0 0 0.00 0.00 0.50 0.50 98 0 98 22 0 1 0 0 0 23 121 0 0 0
May Ruckle 30 H 3 1.50 26.50 0.0 28.00 0.250 0 0.00 1.00 0.00 2.50 0.75 0 0.50 0 0 0 0 0 0 0.00 0.75 4.75 4.75 66 0 66 32 0 1 0 0 0 33 99 0 0 0
May Ruckle 30 H 4 4.50 18.75 18.5 41.75 0.025 0 0.00 0.00 0.25 3.00 7.25 0 0.00 0 0 0 0 0 0 0.00 7.25 10.50 10.50 58 0 58 13 0 0 0 0 1 14 72 0 0 0
May Ruckle 30 H 5 1.50 1.25 0.0 2.75 0.500 0 0.00 0.25 0.00 0.75 71.50 0 0.00 0 0 0 0 0 0 0.00 71.50 72.50 72.50 0 0 0 15 0 0 0 0 1 16 16 5 0 0
May Ruckle 30 H 6 8.25 11.00 0.0 19.25 0.025 0 0.25 0.00 0.00 0.25 15.00 0 0.25 0 0 0 0 0 0 0.25 15.25 15.50 15.50 6 0 6 16 0 1 0 0 0 17 23 9 0 0

Models fit to the data:
1. Zero-inflated Poisson model with a single zero-inflation parameter by applying to all observations (ziformula~1)
2. Zero-inflated Negative Binomial model with NB2 parameterization (variance = \(\mu\)(1 + \(\mu\)/\(\kappa\)), increasing quadratically with the mean)
3. Zero-inflated Negative Binomial model with NB1 parameterization (variance = \(\mu\) + \(\phi\)\(\mu\), increasing linearly with the mean)
4. Poisson model
5. Tweedie model
6. Gaussian model
7. Negative Binomial with NB2 parameterization
8. Negative Binomial with NB1 parameterization

Questions
Do I need to include salinity regime twice in the model syntax? i.e. salinity_regime + (1|salinity_regime/site) YES

Gastropods

Limpets

Because there are so many ‘unknown limpets’ I’m lumping all species together.
Abundance of different limpet species in survey dataset

Abundance of different limpet species in survey dataset

Fitting the models, and using AIC to select best model for the data:

limpet_zipoisson <- glmmTMB(all_limpets_no ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))

limpet_ziNB <- update(limpet_zipoisson, family = nbinom2)

limpet_ziNB1 <- update(limpet_zipoisson, family = nbinom1) # best model 

limpet_poisson <-  glmmTMB(all_limpets_no ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))

limpet_tweedie <- glmmTMB(all_limpets_no ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))

limpet_gaussian <- glmmTMB((all_limpets_no + 1) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))

limpet_NB <- update(limpet_poisson, family = nbinom2)

limpet_NB1 <- update(limpet_poisson, family = nbinom1)
Model Comparison
dAIC df
limpet_ziNB1 0.00000 12
limpet_ziNB 4.31957 12
limpet_NB 5.76179 11
limpet_NB1 12.68349 11
limpet_tweedie 16.21041 12
limpet_gaussian 420.95556 11
limpet_zipoisson 1971.27929 11
limpet_poisson 2357.48378 10

Adding a dispersion formula for salinity regime to account for heterogeneity

limpet_ziNB1_2 <- glmmTMB(all_limpets_no ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, ziformula = ~1, dispformula = ~salinity_regime, family = nbinom1)

#checking assumptions 
limpet_sim_res_2 <- simulateResiduals(limpet_ziNB1_2)
plotResiduals(limpet_sim_res_2, form = survey_data$salinity_regime) # good

Limpet model output: Salinity regime and the interaction with month are significant.
coefficient estimates for model fitted to limpet survey data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.3297650 0.2711593 8.5918671 0.0000000
monthJune -0.1177717 0.2678733 -0.4396546 0.6601873
monthJuly -0.5723872 0.3116878 -1.8364117 0.0662968
monthSeptember -1.7281043 0.3864513 -4.4717257 0.0000078
salinity_regimeH 1.2788265 0.3512972 3.6402978 0.0002723
monthJune:salinity_regimeH 0.2112509 0.3231194 0.6537858 0.5132498
monthJuly:salinity_regimeH 0.5471796 0.3623420 1.5101193 0.1310130
monthSeptember:salinity_regimeH 2.1470353 0.4238716 5.0652966 0.0000004
p values for model fitted to limpet survey data
term statistic df p.value
(Intercept) 73.82018 1 0.0000000
month 22.21336 3 0.0000589
salinity_regime 13.25177 1 0.0002723
month:salinity_regime 27.51028 3 0.0000046
tukey adjusted p values for pairwise comparisons from model fitted to limpet survey data
contrast salinity_regime estimate SE df t.ratio p.value
May - June L 0.1177717 0.2678733 179 0.4396546 0.9715086
May - July L 0.5723872 0.3116878 179 1.8364117 0.2597615
May - September L 1.7281043 0.3864513 179 4.4717257 0.0000806
June - July L 0.4546154 0.3156435 179 1.4402813 0.4760002
June - September L 1.6103326 0.3907160 179 4.1214907 0.0003330
July - September L 1.1557171 0.4194462 179 2.7553405 0.0325067
May - June H -0.0934792 0.1807377 179 -0.5172091 0.9548836
May - July H 0.0252075 0.1853411 179 0.1360061 0.9990990
May - September H -0.4189310 0.1741913 179 -2.4050059 0.0797715
June - July H 0.1186867 0.1823445 179 0.6508927 0.9151248
June - September H -0.3254519 0.1706535 179 -1.9070912 0.2287283
July - September H -0.4441386 0.1747989 179 -2.5408543 0.0571591

limpet abundance from transects in high and low salinity sites over the summer

limpet abundance from transects in high and low salinity sites over the summer

Littorines

littorina spp, sitkana, and littorina total
I’m using total littorines since the maximum number of sitkana is much lower and it’s mostly just Horseshoe Bay that has sitkana. 

littorine_zipoisson <- glmmTMB(littorina_total_no ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))

littorine_ziNB <- update(littorine_zipoisson, family = nbinom2)

littorine_ziNB1 <- update(littorine_zipoisson, family = nbinom1)

littorine_poisson <-  glmmTMB(littorina_total_no ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))

littorine_gaussian <-  glmmTMB((littorina_total_no+1) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))

littorine_tweedie <- update(littorine_poisson, family = tweedie)

littorine_NB <- update(littorine_poisson, family = nbinom2)

littorine_NB1 <- update(littorine_poisson, family = nbinom1)

Littorine model selection using AIC. The best model is the one using a negative binomial distribution, without modelling the zeroes. I needed to add a dispersion formula to model the variance in each salinity region separately.

##                     dAIC    df
## littorine_NB1           0.0 11
## littorine_ziNB1         2.0 12
## littorine_NB           22.7 11
## littorine_ziNB         23.5 12
## littorine_tweedie      67.7 12
## littorine_gaussian    872.5 11
## littorine_zipoisson 11274.6 11
## littorine_poisson   14479.6 10

Adding a dispersion formula to account for heterogeneity in the high and low salinity regions.

# adding in a dispersion formula 
littorine_NB1_2 <- glmmTMB(littorina_total_no ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, dispformula = ~salinity_regime*month, family = nbinom1)

#checking assumptions 
set.seed(10)
littorine_sim_res_2 <- simulateResiduals(littorine_NB1_2)

par(mfrow = c(1,2))
plotResiduals(littorine_sim_res_2, form = survey_data$salinity_regime) # good
plotResiduals(littorine_sim_res_2, form = survey_data$month) # good

Littorine model output:
coefficient estimates for model fitted to littorine survey data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9557704 0.5896000 6.7092439 0.0000000
monthJune -0.2261032 0.6831038 -0.3309939 0.7406491
monthJuly -2.8960237 0.6007762 -4.8204703 0.0000014
monthSeptember -1.0691721 0.7068567 -1.5125725 0.1303883
salinity_regimeH 0.0223592 0.7374873 0.0303181 0.9758133
monthJune:salinity_regimeH -0.5203907 0.7339034 -0.7090725 0.4782795
monthJuly:salinity_regimeH 2.8659941 0.6523130 4.3935874 0.0000111
monthSeptember:salinity_regimeH 0.9981662 0.7608822 1.3118537 0.1895695
p values for model fitted to littorine survey data
term statistic df p.value
month 9.902932 3 0.0194095
salinity_regime 2.380059 1 0.1228929
month:salinity_regime 27.946482 3 0.0000037
tukey adjusted p values for pairwise comparisons from model fitted to littorine survey data
contrast salinity_regime estimate SE df t.ratio p.value
May - June L 0.2261032 0.6831038 174 0.3309939 0.9874522
May - July L 2.8960237 0.6007762 174 4.8204703 0.0000184
May - September L 1.0691721 0.7068567 174 1.5125725 0.4321789
June - July L 2.6699205 0.6935247 174 3.8497843 0.0009466
June - September L 0.8430689 0.7872045 174 1.0709655 0.7076616
July - September L -1.8268516 0.7169325 174 -2.5481502 0.0562093
May - June H 0.7464939 0.2682973 174 2.7823381 0.0302530
May - July H 0.0300296 0.2541264 174 0.1181679 0.9994080
May - September H 0.0710059 0.2815942 174 0.2521569 0.9943548
June - July H -0.7164643 0.2498075 174 -2.8680663 0.0237927
June - September H -0.6754880 0.2777027 174 -2.4324139 0.0748013
July - September H 0.0409763 0.2640371 174 0.1551915 0.9986641

Littorine abundance from transect data in high and low salinity sites over the summer

Littorine abundance from transect data in high and low salinity sites over the summer

Barnacles

species: Balanus, Chthalamus, recruits and total

Balanus models

For percent cover data, I’m rounding to nearest integer.
Best model is negative binomial distribution

balanus_poisson <- glmmTMB(round(balanus_pt) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))

balanus_NB <- update(balanus_poisson, family = nbinom2) 

balanus_NB1 <- update(balanus_poisson, family = nbinom1) #best model 

balanus_gaussian <- glmmTMB((balanus_pt+1) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))

balanus_tweedie <- glmmTMB(balanus_pt ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))

# models with zero-inflation
balanus_zipoisson <- glmmTMB(round(balanus_pt) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))

balanus_ziNB <- update(balanus_zipoisson, family = nbinom2) 

balanus_ziNB1 <- update(balanus_zipoisson, family = nbinom1) 

AICtab(balanus_poisson, balanus_NB, balanus_NB1, balanus_gaussian, balanus_tweedie, balanus_zipoisson, balanus_ziNB, balanus_ziNB1)
##                   dAIC   df
## balanus_NB1          0.0 11
## balanus_ziNB1        2.0 12
## balanus_NB          35.4 11
## balanus_ziNB        37.4 12
## balanus_tweedie     44.5 12
## balanus_gaussian   170.4 11
## balanus_zipoisson 1554.2 11
## balanus_poisson   1817.0 10

Checking model assumptions:

#checking assumptions 
set.seed(6789)
balanus_sim_res <- simulateResiduals(balanus_NB1)

plot(balanus_sim_res)

# plot residuals against time

par(mfrow = c(1,2))
plotResiduals(balanus_sim_res, form = survey_data$month) # different levels of variance in each month 
#plot residuals against salinity regime
plotResiduals(balanus_sim_res, form = survey_data$salinity_regime) # more variance in the low salinity sites

Adding a dispersion formula for salinity regime and month 

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 2.2848, p-value = 0.6995
## alternative hypothesis: true autocorrelation is not 0

Balanus model output:
coefficient estimates for model fitted to Balanus survey data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.9913985 0.3883563 5.1277619 0.0000003
monthJune 0.4464807 0.3207715 1.3918962 0.1639538
monthJuly 0.8889510 0.2854156 3.1145847 0.0018420
monthSeptember 0.7842295 0.3862107 2.0305743 0.0422982
salinity_regimeH 0.3468004 0.5304337 0.6538054 0.5132372
monthJune:salinity_regimeH -0.0164852 0.4184566 -0.0393952 0.9685753
monthJuly:salinity_regimeH 0.0815283 0.3724967 0.2188700 0.8267513
monthSeptember:salinity_regimeH 0.2481843 0.4724092 0.5253588 0.5993338
p values for model fitted to Balanus survey data
term statistic df p.value
month 32.4966405 3 0.0000004
salinity_regime 0.8598830 1 0.3537714
month:salinity_regime 0.4067149 3 0.9388519
tukey adjusted p values for pairwise comparisons from model fitted to Balanus survey data
contrast salinity_regime estimate SE df t.ratio p.value
May - June L -0.4464807 0.3207715 174 -1.3918962 0.5061165
May - July L -0.8889510 0.2854156 174 -3.1145847 0.0114717
May - September L -0.7842295 0.3862107 174 -2.0305743 0.1807860
June - July L -0.4424703 0.2654686 174 -1.6667520 0.3444467
June - September L -0.3377489 0.3717126 174 -0.9086291 0.8002707
July - September L 0.1047215 0.3416693 174 0.3064994 0.9899779
May - June H -0.4299955 0.2687221 174 -1.6001492 0.3812266
May - July H -0.9704793 0.2393569 174 -4.0545283 0.0004367
May - September H -1.0324139 0.2720511 174 -3.7949269 0.0011582
June - July H -0.5404839 0.2128022 174 -2.5398419 0.0574004
June - September H -0.6024184 0.2490077 174 -2.4192765 0.0772000
July - September H -0.0619345 0.2169907 174 -0.2854248 0.9918680

Balanus abundance (% cover) in from transects in high and low salinity sites over the summer

Balanus abundance (% cover) in from transects in high and low salinity sites over the summer

Chthamalus models

Only 1 plot at Rope Site has Chthamalus in it, and no plots at Lions Bay. This is causing convergence issues

chthamalus_poisson <- glmmTMB(round(chthamalus_pt) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = poisson(link = "log")) #convergence problems

chthamalus_NB <- update(chthamalus_poisson, family = nbinom2) 

chthamalus_NB1 <- update(chthamalus_poisson, family = nbinom1) # convergence issues: This often occurs when a model is overparameterized (i.e. the data does not contain information to estimate the parameters).

chthamalus_gaussian <- glmmTMB((chthamalus_pt+1) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))

chthamalus_tweedie <- glmmTMB(chthamalus_pt ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log")) #convergence issues 

# models with zero-inflation
chthamalus_zipoisson <- glmmTMB(round(chthamalus_pt) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log")) #convergence issues

chthamalus_ziNB <- glmmTMB(round(chthamalus_pt) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = nbinom2(link = "log")) # best model 

chthamalus_ziNB1 <- update(chthamalus_ziNB, family = nbinom1) #convergence issues

AICtab(chthamalus_NB, chthamalus_gaussian, chthamalus_ziNB)
##                     dAIC  df
## chthamalus_ziNB       0.0 12
## chthamalus_NB         1.3 11
## chthamalus_gaussian 563.6 11

Checking model assumptions:

Adding a dispersion formula for salinity regime.

Chthamalus model output:
coefficient estimates for model fitted to chthamalus survey data
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.8696123 0.7872713 -1.1045904 0.2693371
monthJune -19.5969277 4449.3512293 -0.0044044 0.9964858
monthJuly -0.6649044 0.9308482 -0.7142995 0.4750420
monthSeptember -0.4624549 0.9193908 -0.5030014 0.6149633
salinity_regimeH 3.8193920 0.9228073 4.1388835 0.0000349
monthJune:salinity_regimeH 19.9082169 4449.3512355 0.0044744 0.9964300
monthJuly:salinity_regimeH 0.5627986 0.9603736 0.5860205 0.5578617
monthSeptember:salinity_regimeH 0.2304269 0.9509023 0.2423245 0.8085287
p values for model fitted to chthamalus survey data
term statistic df p.value
month 6.1112070 3 0.1063233
salinity_regime 22.7604991 1 0.0000018
month:salinity_regime 0.3606038 3 0.9482550
tukey adjusted p values for pairwise comparisons from model fitted to chthamalus survey data
contrast salinity_regime estimate SE df t.ratio p.value
May - June L 19.5969277 4449.3512293 179 0.0044044 1.0000000
May - July L 0.6649044 0.9308482 179 0.7142995 0.8913227
May - September L 0.4624549 0.9193908 179 0.5030014 0.9582756
June - July L -18.9320233 4449.3512080 179 -0.0042550 1.0000000
June - September L -19.1344729 4449.3512047 179 -0.0043005 1.0000000
July - September L -0.2024495 0.8285204 179 -0.2443507 0.9948560
May - June H -0.3112892 0.2332882 179 -1.3343543 0.5423276
May - July H 0.1021057 0.2352839 179 0.4339682 0.9725488
May - September H 0.2320279 0.2418300 179 0.9594674 0.7725801
June - July H 0.4133949 0.2302449 179 1.7954576 0.2789042
June - September H 0.5433171 0.2350804 179 2.3111970 0.0992972
July - September H 0.1299222 0.2369264 179 0.5483653 0.9468890

Chthamalus abundance from transect data in low and high salinity sites over the summer

Chthamalus abundance from transect data in low and high salinity sites over the summer

Invertebrate Plots

Algae

Taxonomic Groups

Reds

Species:
Pyropia, Hildenbrandia, Mastocarpus, Mastocarpus crust, Endocladia, Microcladia, and Polysiphonia
Notes
Mostly Mastocarpus and Pyropia
Not very much red in the Sept/low plots.
Looks like there isn’t a lot at Rope Site.

red_poisson <- glmmTMB(round(reds) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))

red_NB <- update(red_poisson, family = nbinom2) 

red_NB1 <- update(red_poisson, family = nbinom1) #best model 

red_gaussian <- glmmTMB((reds+1) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))

red_tweedie <- glmmTMB(reds ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))

# models with zero-inflation
red_zipoisson <- glmmTMB(round(reds) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))

red_ziNB <- update(red_zipoisson, family = nbinom2) #convergence issues

red_ziNB1 <- update(red_zipoisson, family = nbinom1)

There is basically no red algae at Rope Site, which is making it very difficult to fit a model that meets the assumption of homogeneity of variance.

##               dAIC  df
## red_NB1         0.0 11
## red_ziNB1       2.0 12
## red_NB         19.6 11
## red_tweedie    78.9 12
## red_zipoisson 508.9 11
## red_gaussian  509.1 11
## red_poisson   792.1 10

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 1.1046, p-value = 0.2805
## alternative hypothesis: true autocorrelation is not 0

Adding a dispersion parameter for salinity_regime:

Adding dispersion parameter for site instead:

Red Algae model output:
coefficient estimates for model fitted to red algae survey data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.3955417 0.3310490 4.2155134 0.0000249
monthJune -0.5682088 0.4505585 -1.2611212 0.2072652
monthJuly -0.3527778 0.4088760 -0.8627989 0.3882481
monthSeptember -0.9547281 0.5077820 -1.8801931 0.0600818
salinity_regimeH 0.6056305 0.3783527 1.6007037 0.1094426
monthJune:salinity_regimeH 0.7737543 0.4948812 1.5635154 0.1179314
monthJuly:salinity_regimeH 0.4175858 0.4610854 0.9056582 0.3651168
monthSeptember:salinity_regimeH 0.6503174 0.5589100 1.1635457 0.2446082
p values for model fitted to red algae survey data
term statistic df p.value
month 6.490347 3 0.0900440
salinity_regime 12.655822 1 0.0003744
month:salinity_regime 2.850151 3 0.4153106
tukey adjusted p values for pairwise comparisons from model fitted to red algae survey data
contrast salinity_regime estimate SE df t.ratio p.value
May - June L 0.5682088 0.4505585 176 1.2611212 0.5888808
May - July L 0.3527778 0.4088760 176 0.8627989 0.8239858
May - September L 0.9547281 0.5077820 176 1.8801931 0.2402994
June - July L -0.2154311 0.4594790 176 -0.4688595 0.9657853
June - September L 0.3865193 0.5523964 176 0.6997136 0.8970712
July - September L 0.6019503 0.5177022 176 1.1627347 0.6511123
May - June H -0.2055455 0.2047147 176 -1.0040583 0.7472381
May - July H -0.0648080 0.2131860 176 -0.3039975 0.9902161
May - September H 0.3044107 0.2335373 176 1.3034778 0.5619437
June - July H 0.1407375 0.2003558 176 0.7024377 0.8960095
June - September H 0.5099562 0.2225810 176 2.2911036 0.1040051
July - September H 0.3692187 0.2288667 176 1.6132477 0.3738178

Red Algae (percent cover) from transect data in low and high salinity sites over the summer

Red Algae (percent cover) from transect data in low and high salinity sites over the summer

Greens

Species: ulva_pt, urospora_pt, acrosiphonia_pt

Notes
Mostly Ulva 

green_poisson <- glmmTMB(round(greens) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))

green_NB <- update(green_poisson, family = nbinom2) 

green_NB1 <- update(green_poisson, family = nbinom1) #best model 

green_gaussian <- glmmTMB((greens+1) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))

green_tweedie <- glmmTMB(greens ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))

# models with zero-inflation
green_zipoisson <- glmmTMB(round(greens) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))

green_ziNB <- update(green_zipoisson, family = nbinom2) 

green_ziNB1 <- update(green_zipoisson, family = nbinom1) #convergence issues

Checking model assumptions:

##                 dAIC  df
## green_NB1         0.0 11
## green_NB         11.3 11
## green_ziNB       13.3 12
## green_tweedie    88.0 12
## green_zipoisson 462.5 11
## green_gaussian  610.4 11
## green_poisson   951.5 10

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 2.1239, p-value = 0.8461
## alternative hypothesis: true autocorrelation is not 0

Green Algae model summary:
coefficient estimates for model fitted to green algae survey data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.3910871 0.5332449 0.7334101 0.4633084
monthJune 0.2256104 0.4656709 0.4844845 0.6280420
monthJuly 0.8711570 0.4315772 2.0185428 0.0435348
monthSeptember 1.9120327 0.4043576 4.7285685 0.0000023
salinity_regimeH 0.3165282 0.7172861 0.4412858 0.6590061
monthJune:salinity_regimeH -0.2971707 0.6931032 -0.4287539 0.6681023
monthJuly:salinity_regimeH -0.7566436 0.6465476 -1.1702829 0.2418871
monthSeptember:salinity_regimeH -0.9512348 0.5879712 -1.6178253 0.1057002
p values for model fitted to green algae survey data
term statistic df p.value
month 39.3806196 3 0.0000000
salinity_regime 0.3592205 1 0.5489395
month:salinity_regime 3.2021357 3 0.3614974
tukey adjusted p values for pairwise comparisons from model fitted to green algae survey data
contrast salinity_regime estimate SE df t.ratio p.value
May - June L -0.2256104 0.4656709 181 -0.4844845 0.9624619
May - July L -0.8711570 0.4315772 181 -2.0185428 0.1849684
May - September L -1.9120327 0.4043576 181 -4.7285685 0.0000268
June - July L -0.6455467 0.4016562 181 -1.6072121 0.3771338
June - September L -1.6864224 0.3727681 181 -4.5240523 0.0000644
July - September L -1.0408757 0.3202564 181 -3.2501328 0.0074540
May - June H 0.0715604 0.5133258 181 0.1394053 0.9990301
May - July H -0.1145134 0.4815480 181 -0.2378026 0.9952531
May - September H -0.9607980 0.4320781 181 -2.2236674 0.1207383
June - July H -0.1860737 0.4980846 181 -0.3735786 0.9821648
June - September H -1.0323583 0.4450468 181 -2.3196621 0.0973494
July - September H -0.8462846 0.4105998 181 -2.0610932 0.1699512

Green Algae from transects in low and high salinity sites over the summer

Green Algae from transects in low and high salinity sites over the summer

Browns

Species: fucus, melanosiphon, leathesia
Notes
Mostly fucus

brown_poisson <- glmmTMB(round(browns) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))

brown_NB <- update(brown_poisson, family = nbinom2) 

brown_NB1 <- update(brown_poisson, family = nbinom1) #best model 

brown_gaussian <- glmmTMB((browns+1) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))

brown_tweedie <- glmmTMB(browns ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))

# models with zero-inflation
brown_zipoisson <- glmmTMB(round(browns) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))

brown_ziNB <- update(brown_zipoisson, family = nbinom2) 

brown_ziNB1 <- update(brown_zipoisson, family = nbinom1)

Checking assumptions:

##                 dAIC   df
## brown_NB1          0.0 11
## brown_ziNB1        1.5 12
## brown_NB          22.4 11
## brown_ziNB        24.4 12
## brown_tweedie     51.7 12
## brown_gaussian   348.9 11
## brown_zipoisson 2856.6 11
## brown_poisson   3237.8 10

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 2.1239, p-value = 0.8461
## alternative hypothesis: true autocorrelation is not 0

Adding a dispersion parameter for salinity_regime:

Adding dispersion parameter for site instead:

Brown algae model summary:
coefficient estimates for model fitted to littorine survey data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.0445515 0.1928951 15.7834597 0.0000000
monthJune 0.1421746 0.2523742 0.5633484 0.5731977
monthJuly 0.6544828 0.2303375 2.8414082 0.0044915
monthSeptember 0.7550187 0.2290348 3.2965238 0.0009789
salinity_regimeH 0.2522594 0.3170038 0.7957612 0.4261708
monthJune:salinity_regimeH -0.4583030 0.4169180 -1.0992641 0.2716529
monthJuly:salinity_regimeH -0.6180043 0.3877060 -1.5940025 0.1109355
monthSeptember:salinity_regimeH -0.6116013 0.3798392 -1.6101584 0.1073633
p values for model fitted to brown algae survey data
term statistic df p.value
month 15.399269 3 0.0015054
salinity_regime 1.207230 1 0.2718814
month:salinity_regime 3.324017 3 0.3443139
tukey adjusted p values for pairwise comparisons from model fitted to brown algae survey data
contrast salinity_regime estimate SE df t.ratio p.value
May - June L -0.1421746 0.2523742 176 -0.5633484 0.9427681
May - July L -0.6544828 0.2303375 176 -2.8414082 0.0256300
May - September L -0.7550187 0.2290348 176 -3.2965238 0.0064487
June - July L -0.5123083 0.2153927 176 -2.3784844 0.0850085
June - September L -0.6128441 0.2305027 176 -2.6587281 0.0422054
July - September L -0.1005358 0.2071795 176 -0.4852596 0.9622891
May - June H 0.3161285 0.3318553 176 0.9526094 0.7763957
May - July H -0.0364786 0.3118662 176 -0.1169686 0.9994258
May - September H -0.1434174 0.3030196 176 -0.4732940 0.9648602
June - July H -0.3526070 0.3302644 176 -1.0676508 0.7096561
June - September H -0.4595458 0.3224557 176 -1.4251441 0.4853869
July - September H -0.1069388 0.3018381 176 -0.3542919 0.9847048

Taxonomic Algae Plot

#### Functional Algal Groupings Grazed (sheets and filaments) vs. Non grazed (leathery, calcareous, crustose). Is there a better term for this?

Grazed Filaments, sheets Ulva, Pyropia, Polysiphonia, Diatoms, Urospora, Acrosiphonia

Non-grazed Turfs, crusts Fucus, Mastocarpus, Mastocarpus crust, Hildenbrandia, Endocladia, Melanosiphon (not sure if this is the right grouping), leathesia Endocladia and Mastocarpus both part of Lottia diet though?

Grazed Algae
grazed_algae_poisson <- glmmTMB(round(grazed) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))

grazed_algae_NB <- update(grazed_algae_poisson, family = nbinom2) # Best model 

grazed_algae_NB1 <- update(grazed_algae_poisson, family = nbinom1) 

grazed_algae_gaussian <- glmmTMB((grazed+1) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))

grazed_algae_tweedie <- glmmTMB(grazed ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))

# models with zero-inflation
grazed_algae_zipoisson <- glmmTMB(round(grazed) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))

grazed_algae_ziNB <- update(grazed_algae_zipoisson, family = nbinom2)

grazed_algae_ziNB1 <- update(grazed_algae_zipoisson, family = nbinom1)

Grazed algae: checking model assumptions

##                        dAIC   df
## grazed_algae_NB           0.0 11
## grazed_algae_NB1          1.4 11
## grazed_algae_ziNB         2.0 12
## grazed_algae_ziNB1        2.4 12
## grazed_algae_tweedie     65.0 12
## grazed_algae_gaussian   576.2 11
## grazed_algae_zipoisson  781.4 11
## grazed_algae_poisson   1330.9 10

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 2.734, p-value = 0.3787
## alternative hypothesis: true autocorrelation is not 0

Grazed algae model summary
coefficient estimates for model fitted to grazed algae survey data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.2750823 0.6570517 0.4186615 0.6754635
monthJune 0.3417898 0.5474157 0.6243698 0.5323847
monthJuly 1.0713042 0.5411664 1.9796208 0.0477462
monthSeptember 2.0280985 0.5326420 3.8076204 0.0001403
salinity_regimeH 1.4646335 0.9147464 1.6011361 0.1093468
monthJune:salinity_regimeH -0.7042222 0.7381904 -0.9539844 0.3400915
monthJuly:salinity_regimeH -1.6986907 0.7425933 -2.2875116 0.0221660
monthSeptember:salinity_regimeH -2.0856043 0.7412548 -2.8136132 0.0048988
p values for model fitted to grazed algae survey data
term statistic df p.value
month 10.2155670 3 0.0168199
salinity_regime 0.1438865 1 0.7044475
month:salinity_regime 9.8419923 3 0.0199581
tukey adjusted p values for pairwise comparisons from model fitted to grazed algae survey data
contrast salinity_regime estimate SE df t.ratio p.value
May - June L -0.3417898 0.5474157 181 -0.6243698 0.9241470
May - July L -1.0713042 0.5411664 181 -1.9796208 0.1995108
May - September L -2.0280985 0.5326420 181 -3.8076204 0.0010930
June - July L -0.7295144 0.5527433 181 -1.3198070 0.5515370
June - September L -1.6863087 0.5079056 181 -3.3201222 0.0059442
July - September L -0.9567942 0.5318073 181 -1.7991371 0.2771120
May - June H 0.3624323 0.4956890 181 0.7311688 0.8844714
May - July H 0.6273865 0.5088229 181 1.2330153 0.6067162
May - September H 0.0575059 0.5131515 181 0.1120641 0.9994949
June - July H 0.2649541 0.5220343 181 0.5075416 0.9572101
June - September H -0.3049265 0.5198372 181 -0.5865807 0.9360315
July - September H -0.5698806 0.5228355 181 -1.0899806 0.6961113

Non grazed algae
non_grazed_algae_poisson <- glmmTMB(round(non_grazed) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))

non_grazed_algae_NB <- update(non_grazed_algae_poisson, family = nbinom2) # Best model 

non_grazed_algae_NB1 <- update(non_grazed_algae_poisson, family = nbinom1) 

non_grazed_algae_gaussian <- glmmTMB((non_grazed+1) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))

non_grazed_algae_tweedie <- glmmTMB(non_grazed ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))

# models with zero-inflation
non_grazed_algae_zipoisson <- glmmTMB(round(non_grazed) ~ month * salinity_regime +  (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))

non_grazed_algae_ziNB <- update(non_grazed_algae_zipoisson, family = nbinom2)
# convergence issues 

non_grazed_algae_ziNB1 <- update(non_grazed_algae_zipoisson, family = nbinom1)

AICtab(non_grazed_algae_poisson, non_grazed_algae_NB, non_grazed_algae_NB1, non_grazed_algae_gaussian, non_grazed_algae_tweedie, non_grazed_algae_zipoisson, non_grazed_algae_ziNB1)
##                            dAIC   df
## non_grazed_algae_NB1          0.0 11
## non_grazed_algae_ziNB1        1.5 12
## non_grazed_algae_NB          22.2 11
## non_grazed_algae_tweedie     48.3 12
## non_grazed_algae_gaussian   349.2 11
## non_grazed_algae_zipoisson 2852.3 11
## non_grazed_algae_poisson   3230.9 10

Checking model assumptions:

Adding a dispersion parameter for salinity regime:

Adding a dispersion formula for site instead: 

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 1.8743, p-value = 0.8443
## alternative hypothesis: true autocorrelation is not 0

Non-Grazed Algae model summary:
coefficient estimates for model fitted to Non-grazed algae survey data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.0289719 0.1924510 15.7389283 0.0000000
monthJune 0.1414859 0.2524612 0.5604265 0.5751886
monthJuly 0.6639990 0.2299227 2.8879226 0.0038780
monthSeptember 0.7695971 0.2285388 3.3674681 0.0007586
salinity_regimeH 0.2678394 0.3167337 0.8456296 0.3977594
monthJune:salinity_regimeH -0.4576148 0.4169706 -1.0974750 0.2724338
monthJuly:salinity_regimeH -0.6275224 0.3874597 -1.6195811 0.1053223
monthSeptember:salinity_regimeH -0.6261807 0.3795403 -1.6498400 0.0989757
p values for model fitted to Non-grazed algae survey data
term statistic df p.value
month 15.958681 3 0.0011563
salinity_regime 1.111575 1 0.2917398
month:salinity_regime 3.459912 3 0.3259987
tukey adjusted p values for pairwise comparisons from model fitted to Non-grazed algae survey data
contrast salinity_regime estimate SE df t.ratio p.value
May - June L -0.1414859 0.2524612 176 -0.5604265 0.9435853
May - July L -0.6639990 0.2299227 176 -2.8879226 0.0224567
May - September L -0.7695971 0.2285388 176 -3.3674681 0.0051136
June - July L -0.5225130 0.2156012 176 -2.4235164 0.0763768
June - September L -0.6281111 0.2305096 176 -2.7248811 0.0353643
July - September L -0.1055981 0.2066763 176 -0.5109348 0.9563991
May - June H 0.3161289 0.3318551 176 0.9526111 0.7763948
May - July H -0.0364766 0.3118663 176 -0.1169623 0.9994259
May - September H -0.1434164 0.3030195 176 -0.4732909 0.9648608
June - July H -0.3526055 0.3302644 176 -1.0676460 0.7096591
June - September H -0.4595452 0.3224556 176 -1.4251425 0.4853879
July - September H -0.1069398 0.3018382 176 -0.3542950 0.9847044

Abundance of non-grazed algae from transects in low and high salinity sites over the summer

Abundance of non-grazed algae from transects in low and high salinity sites over the summer

Survey functional algae plot and anova table

## Field Experiment Data Notes
Only using the control and exclusion plots.
Not modelling limpets b/c they will be different due to treatment.
Do we just assume the herbivores that were removed from the exclusion plots were not having an impact on the community? i.e. they weren’t there for long enough?
Waiting to hear back from Theraesa about:
Additional Limpets in ring (Removed) two of the control plots show 34 and 1 paradigitalis in this column but 0s under the paradigitalis column.
Total Limpets (Non-Pelta removed from rings) some of the control plots have other limpet species as removed? Littorina Removed from Ring (#) some of the control plots have littorines removed?

Species: Balanus, Chthamalus, limpets, littorina, mytilus, amphipods diatoms, pyropia, fucus, ulva,

Limpets

Check and see whether the exclusion plots had less limpets than the controls. We expect the interaction to be significant as well. I’m using data from May, June and July for this as previous months will have had an impact on July’s communities as well.

##                       dAIC   df
## exp_limpets_NB           0.0 7 
## exp_limpets_NB1         10.4 7 
## exp_limpets_tweedie     29.4 8 
## exp_limpets_zipoisson   81.0 7 
## exp_limpets_poisson    391.0 6 
## exp_limpets_gaussian  1079.8 7

Limpets model output:
coefficient estimates for model fitted to Limpet field experiment data for May, June and July
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.2938817 0.9720340 -3.3886488 0.0007024
treatmentE 0.6051422 0.9769437 0.6194238 0.5356372
regionHigh 3.2174139 1.2048968 2.6702817 0.0075788
treatmentE:regionHigh -1.8579827 1.1802315 -1.5742529 0.1154290
p values for model fitted to Limpet field experiment data
term statistic df p.value
treatment 1.519823 1 0.2176463
region 4.699369 1 0.0301737
treatment:region 2.478272 1 0.1154290
tukey adjusted p values for pairwise comparisons from model fitted to Limpet field experiment data
contrast region estimate SE df t.ratio p.value
C - E Low -0.6051422 0.9769437 245 -0.6194238 0.5362128
C - E High 1.2528406 0.6580842 245 1.9037694 0.0581120

Balanus

Question How were barnacles counted? How can we have non-integers? 

exp_balanus_poisson <- glmmTMB(round(balanus_no) ~ treatment * region +  (1|region/site), data = exp_data, family = poisson(link = "log"))

exp_balanus_NB <- update(exp_balanus_poisson, family = nbinom2) 

exp_balanus_NB1 <- update(exp_balanus_poisson, family = nbinom1) 

exp_balanus_gaussian <- glmmTMB((balanus_no+1) ~ treatment * region +  (1|region/site), data = exp_data, family = gaussian(link = "log"))

exp_balanus_tweedie <- glmmTMB(balanus_no ~ treatment * region +  (1|region/site), data = exp_data, family = tweedie(link = "log")) # Best model 

# models with zero-inflation
exp_balanus_zipoisson <- glmmTMB(round(balanus_no) ~ treatment * region +  (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))

exp_balanus_ziNB <- update(exp_balanus_zipoisson, family = nbinom2) # convergence issues

exp_balanus_ziNB1 <- update(exp_balanus_zipoisson, family = nbinom1) # convergence issues

Checking model assumptions:

##                       dAIC    df
## exp_balanus_tweedie       0.0 8 
## exp_balanus_NB1           2.4 7 
## exp_balanus_gaussian     10.0 7 
## exp_balanus_NB           20.4 7 
## exp_balanus_poisson   13064.6 6 
## exp_balanus_zipoisson 13066.6 7

Adding a dispersion formula for trmt*region:

Trying to fit dispersion formula to treatment*site but it’s causing convergence issues.

What if I don’t include site as a random effect, but as a fixed effect instead?

Adding a dispersion parameter for treatment*site

Balanus model output:
coefficient estimates for model fitted to Balanus field experiment data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.1572444 0.3390990 15.2086678 0.0000000
treatmentexclusion 0.6170414 0.3796994 1.6250790 0.1041457
siteRope Site 2 2.0630250 0.3410859 6.0484033 0.0000000
siteHorseshoe Bay 0.8385693 0.4697243 1.7852373 0.0742228
siteRuckle 1.7846913 0.4240018 4.2091596 0.0000256
siteEagle Cove 2.3308623 0.3667312 6.3557795 0.0000000
siteHailstorm 2 1.2292992 0.4255567 2.8886850 0.0038686
treatmentexclusion:siteRope Site 2 -0.6117046 0.3827412 -1.5982199 0.1099941
treatmentexclusion:siteHorseshoe Bay -0.0399731 0.5254315 -0.0760767 0.9393581
treatmentexclusion:siteRuckle -0.8247594 0.4626999 -1.7824933 0.0746688
treatmentexclusion:siteEagle Cove -0.8251589 0.4296878 -1.9203683 0.0548114
treatmentexclusion:siteHailstorm 2 -0.2195158 0.4683627 -0.4686877 0.6392928
p values for model fitted to Balanus field experiment data
term statistic df p.value
(Intercept) 231.303576 1 0.0000000
treatment 2.640882 1 0.1041457
site 65.318524 5 0.0000000
treatment:site 8.834319 5 0.1158579
tukey adjusted p values for pairwise comparisons from model fitted to Balanus field experiment data
contrast estimate SE df t.ratio p.value
Lions Bay 2 - Rope Site 2 -1.7571727 0.1913706 59 -9.182041 0.0000000
Lions Bay 2 - Horseshoe Bay -0.8185828 0.2627158 59 -3.115850 0.0321635
Lions Bay 2 - Ruckle -1.3723116 0.2313499 59 -5.931757 0.0000025
Lions Bay 2 - Eagle Cove -1.9182829 0.2148439 59 -8.928728 0.0000000
Lions Bay 2 - Hailstorm 2 -1.1195413 0.2341813 59 -4.780660 0.0001691
Rope Site 2 - Horseshoe Bay 0.9385899 0.1831843 59 5.123747 0.0000495
Rope Site 2 - Ruckle 0.3848610 0.1343864 59 2.863840 0.0612430
Rope Site 2 - Eagle Cove -0.1611102 0.1034158 59 -1.557888 0.6288773
Rope Site 2 - Hailstorm 2 0.6376314 0.1392042 59 4.580547 0.0003403
Horseshoe Bay - Ruckle -0.5537288 0.2246254 59 -2.465121 0.1514418
Horseshoe Bay - Eagle Cove -1.0997001 0.2075854 59 -5.297580 0.0000263
Horseshoe Bay - Hailstorm 2 -0.3009585 0.2275405 59 -1.322659 0.7713864
Ruckle - Eagle Cove -0.5459712 0.1661171 59 -3.286664 0.0202130
Ruckle - Hailstorm 2 0.2527703 0.1904675 59 1.327104 0.7689080
Eagle Cove - Hailstorm 2 0.7987416 0.1700382 59 4.697423 0.0002266

Balanus (percent cover) in control vs. exclusion plots in low and high salinity sites

Balanus (percent cover) in control vs. exclusion plots in low and high salinity sites

Chthamalus

exp_chthamalus_poisson <- glmmTMB(round(chthamalus_no) ~ treatment * region +  (1|region/site), data = exp_data, family = poisson(link = "log"))

exp_chthamalus_NB <- update(exp_chthamalus_poisson, family = nbinom2) 

exp_chthamalus_NB1 <- update(exp_chthamalus_poisson, family = nbinom1) 

exp_chthamalus_gaussian <- glmmTMB((chthamalus_no+1) ~ treatment * region +  (1|region/site), data = exp_data, family = gaussian(link = "log"))

exp_chthamalus_tweedie <- glmmTMB(chthamalus_no ~ treatment * region +  (1|region/site), data = exp_data, family = tweedie(link = "log")) 

# models with zero-inflation
exp_chthamalus_zipoisson <- glmmTMB(round(chthamalus_no) ~ treatment * region +  (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))

exp_chthamalus_ziNB <- update(exp_chthamalus_zipoisson, family = nbinom2) 

exp_chthamalus_ziNB1 <- update(exp_chthamalus_zipoisson, family = nbinom1) # best model

Model Selection and checking assumptions:

##                          dAIC   df
## exp_chthamalus_ziNB1        0.0 8 
## exp_chthamalus_tweedie      3.2 8 
## exp_chthamalus_NB1          5.7 7 
## exp_chthamalus_ziNB        19.0 8 
## exp_chthamalus_NB          27.8 7 
## exp_chthamalus_gaussian   402.8 7 
## exp_chthamalus_zipoisson 2222.5 7 
## exp_chthamalus_poisson   3612.3 6

coefficient estimates for model fitted to Chthamalus field experiment data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.3549717 1.0679766 2.2050780 0.0274486
treatmentexclusion 0.0791074 0.3619432 0.2185629 0.8269905
regionHigh 2.3623231 1.4142626 1.6703567 0.0948488
treatmentexclusion:regionHigh -2.2124492 0.5343787 -4.1402269 0.0000347
p values for model fitted to Chthamalus field experiment data
term statistic df p.value
treatment 12.613247 1 0.0003830
region 1.441596 1 0.2298812
treatment:region 17.141479 1 0.0000347
tukey adjusted p values for pairwise comparisons from model fitted to Chthamalus field experiment data
contrast region estimate SE df t.ratio p.value
control - exclusion Low -0.0791074 0.3619432 76 -0.2185629 0.8275764
control - exclusion High 2.1333418 0.3913361 76 5.4514317 0.0000006

Chthamalus (percent cover) in control vs. exclusion plots in low and high salinity sites

Chthamalus (percent cover) in control vs. exclusion plots in low and high salinity sites

Barnacle plot and anova tables

Algae

## # A tibble: 6 × 2
##   site              n
##   <fct>         <int>
## 1 Lions Bay 2      14
## 2 Rope Site 2      14
## 3 Horseshoe Bay    14
## 4 Ruckle           11
## 5 Eagle Cove        8
## 6 Hailstorm 2      12

Fit models to all algae pooled together:

exp_algae_poisson <- glmmTMB(round(algae) ~ treatment * region +  (1|region/site), data = exp_data, family = poisson(link = "log"))

exp_algae_NB <- update(exp_algae_poisson, family = nbinom2) 

exp_algae_NB1 <- update(exp_algae_poisson, family = nbinom1) 

exp_algae_gaussian <- glmmTMB((algae + 1) ~ treatment * region +  (1|region/site), data = exp_data, family = gaussian(link = "log"))

exp_algae_tweedie <- glmmTMB(algae ~ treatment * region +  (1|region/site), data = exp_data, family = tweedie(link = "log"))

# models with zero-inflation
exp_algae_zipoisson <- glmmTMB(round(algae) ~ treatment * region +  (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))

exp_algae_ziNB <- update(exp_algae_zipoisson, family = nbinom2) 

exp_algae_ziNB1 <- update(exp_algae_zipoisson, family = nbinom1) #best model

Checking model assumptions:

##                     dAIC   df
## exp_algae_ziNB1        0.0 8 
## exp_algae_NB1          5.8 7 
## exp_algae_tweedie     15.8 8 
## exp_algae_ziNB        34.7 8 
## exp_algae_NB          41.9 7 
## exp_algae_gaussian    76.5 7 
## exp_algae_zipoisson  764.4 7 
## exp_algae_poisson   1117.8 6

Adding a dispersion parameter for region:

model output:
coefficient estimates for model fitted to Algae field experiment data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.045601 0.2073427 19.5116589 0.0000000
treatmentexclusion -0.115691 0.1551350 -0.7457439 0.4558221
regionHigh -2.179333 0.4570558 -4.7681986 0.0000019
treatmentexclusion:regionHigh 2.274326 0.4043492 5.6246579 0.0000000
p values for model fitted to Algae field experiment data
term statistic df p.value
treatment 2.3728287 1 0.1234632
region 0.8474256 1 0.3572816
treatment:region 31.6367762 1 0.0000000
tukey adjusted p values for pairwise comparisons from model fitted to Algae field experiment data
contrast region estimate SE df t.ratio p.value
control - exclusion Low 0.115691 0.1551350 75 0.7457439 0.4581531
control - exclusion High -2.158635 0.3731297 75 -5.7852129 0.0000002

Algae in control vs. exclusion plots in low and high salinity sites

Algae in control vs. exclusion plots in low and high salinity sites

Taxonomic groupings

Red Algae

Almost all Mastocarpus crust with a very small amount of Pyropia 

exp_red_poisson <- glmmTMB(round(reds) ~ treatment * region +  (1|region/site), data = exp_data, family = poisson(link = "log"))

exp_red_NB <- update(exp_red_poisson, family = nbinom2) 

exp_red_NB1 <- update(exp_red_poisson, family = nbinom1) 

exp_red_gaussian <- glmmTMB((reds + 1) ~ treatment * region +  (1|region/site), data = exp_data, family = gaussian(link = "log"))

exp_red_tweedie <- glmmTMB(reds ~ treatment * region +  (1|region/site), data = exp_data, family = tweedie(link = "log")) # best model 

# models with zero-inflation
exp_red_zipoisson <- glmmTMB(round(reds) ~ treatment * region +  (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))

exp_red_ziNB <- update(exp_red_zipoisson, family = nbinom2) 

exp_red_ziNB1 <- update(exp_red_zipoisson, family = nbinom1) 
##                   dAIC  df
## exp_red_tweedie     0.0 8 
## exp_red_zipoisson   4.4 7 
## exp_red_ziNB1       4.9 8 
## exp_red_ziNB        6.1 8 
## exp_red_NB1         9.8 7 
## exp_red_NB         16.7 7 
## exp_red_poisson    73.3 6 
## exp_red_gaussian  172.5 7

coefficient estimates for model fitted to Red Algae field experiment data
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.0234567 0.7880334 -0.0297661 0.9762536
treatmentexclusion -0.3085752 0.4686554 -0.6584266 0.5102640
regionHigh -0.4996745 1.0881564 -0.4591937 0.6460951
treatmentexclusion:regionHigh -0.2724478 0.8529393 -0.3194223 0.7494063
p values for model fitted to Red Algae field experiment data
term statistic df p.value
treatment 0.9961414 1 0.3182460
region 0.3475700 1 0.5554919
treatment:region 0.1020306 1 0.7494063
tukey adjusted p values for pairwise comparisons from model fitted to Red Algae field experiment data
contrast region estimate SE df t.ratio p.value
control - exclusion Low 0.3085752 0.4686554 76 0.6584266 0.5122538
control - exclusion High 0.5810230 0.7126650 76 0.8152820 0.4174586

Red Algae in control vs. limpet exclusion plots in high and low salinity sites

Red Algae in control vs. limpet exclusion plots in high and low salinity sites

Green Algae

Only green algae species is Ulva
Only region is significant, because the confidence intervals are very large.
How many plots, and in which sites, have green algae?

## # A tibble: 6 × 2
##   site              n
##   <fct>         <int>
## 1 Lions Bay 2      14
## 2 Rope Site 2      14
## 3 Horseshoe Bay    12
## 4 Ruckle            8
## 5 Eagle Cove        6
## 6 Hailstorm 2       5

exp_green_poisson <- glmmTMB(round(greens) ~ treatment * region +  (1|region/site), data = exp_data, family = poisson(link = "log"))

exp_green_NB <- update(exp_green_poisson, family = nbinom2) 

exp_green_NB1 <- update(exp_green_poisson, family = nbinom1) 

exp_green_gaussian <- glmmTMB((greens+1) ~ treatment * region +  (1|region/site), data = exp_data, family = gaussian(link = "log"))

exp_green_tweedie <- glmmTMB(greens ~ treatment * region +  (1|region/site), data = exp_data, family = tweedie(link = "log")) 

# models with zero-inflation
exp_green_zipoisson <- glmmTMB(round(greens) ~ treatment * region +  (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))

exp_green_ziNB <- update(exp_green_zipoisson, family = nbinom2) 

exp_green_ziNB1 <- update(exp_green_zipoisson, family = nbinom1) # best model 
##                     dAIC   df
## exp_green_ziNB1        0.0 8 
## exp_green_tweedie     11.4 8 
## exp_green_NB1         14.1 7 
## exp_green_ziNB        43.3 8 
## exp_green_NB          79.1 7 
## exp_green_gaussian   166.8 7 
## exp_green_zipoisson  539.5 7 
## exp_green_poisson   1310.0 6

Adding a dispersion formula for region:

Green Algae Model outputs:
coefficient estimates for model fitted to Green Algae field experiment data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.0687806 0.1415476 28.7449540 0.0000000
treatmentexclusion -0.1154469 0.1417591 -0.8143877 0.4154229
regionHigh -3.5833048 0.8293180 -4.3207852 0.0000155
treatmentexclusion:regionHigh 3.5105808 0.8378930 4.1897723 0.0000279
p values for model fitted to Green Algae field experiment data
term statistic df p.value
(Intercept) 826.2723828 1 0.0000000
treatment 0.6632273 1 0.4154229
region 18.6691845 1 0.0000155
treatment:region 17.5541920 1 0.0000279
tukey adjusted p values for pairwise comparisons from model fitted to Green Algae field experiment data
contrast region estimate SE df t.ratio p.value
control - exclusion Low 0.1154469 0.1417591 75 0.8143877 0.4180012
control - exclusion High -3.3951339 0.8248749 75 -4.1159380 0.0000980

Brown Algae

Only brown algae species is Fucus
How many plots, and in which sites, have brown algae?
No brown algae at Ruckle Park

## # A tibble: 5 × 2
##   site              n
##   <fct>         <int>
## 1 Lions Bay 2       3
## 2 Rope Site 2       2
## 3 Horseshoe Bay     3
## 4 Eagle Cove        2
## 5 Hailstorm 2       1

exp_brown_poisson <- glmmTMB(round(browns) ~ treatment * region +  (1|region/site), data = exp_data, family = poisson(link = "log"))

exp_brown_NB <- update(exp_brown_poisson, family = nbinom2) 

exp_brown_NB1 <- update(exp_brown_poisson, family = nbinom1) 

exp_brown_gaussian <- glmmTMB((browns + 1) ~ treatment * region +  (1|region/site), data = exp_data, family = gaussian(link = "log"))

exp_brown_tweedie <- glmmTMB(browns ~ treatment * region +  (1|region/site), data = exp_data, family = tweedie(link = "log")) # best model 

# models with zero-inflation
exp_brown_zipoisson <- glmmTMB(round(browns) ~ treatment * region +  (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))

exp_brown_ziNB <- update(exp_brown_zipoisson, family = nbinom2) #convergence issues

exp_brown_ziNB1 <- update(exp_brown_zipoisson, family = nbinom1)

Checking assumptions:

##                     dAIC df
## exp_brown_tweedie    0.0 8 
## exp_brown_zipoisson  7.3 7 
## exp_brown_ziNB1      9.3 8 
## exp_brown_NB1       10.8 7 
## exp_brown_NB        11.0 7 
## exp_brown_poisson   13.5 6 
## exp_brown_gaussian  25.9 7

Brown Algae Model outputs:
coefficient estimates for model fitted to Brown Algae field experiment data
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.8934872 4.104650e-01 -2.1767684 0.0294979
treatmentexclusion -1.0116009 7.677946e-01 -1.3175411 0.1876573
regionHigh -1.2992830 8.454432e-01 -1.5368069 0.1243406
treatmentexclusion:regionHigh -21.3472632 1.932017e+04 -0.0011049 0.9991184
p values for model fitted to Brown Algae field experiment data
term statistic df p.value
(Intercept) 4.7383208 1 0.0294979
treatment 1.7359146 1 0.1876573
region 2.3617753 1 0.1243406
treatment:region 0.0000012 1 0.9991184
tukey adjusted p values for pairwise comparisons from model fitted to Brown Algae field experiment data
contrast region estimate SE df t.ratio p.value
control - exclusion Low 1.011601 7.677946e-01 76 1.3175411 0.1916155
control - exclusion High 22.358864 1.932017e+04 76 0.0011573 0.9990797

Brown Algae in control vs. limpet exclusion plots in high and low salinity sites

Brown Algae in control vs. limpet exclusion plots in high and low salinity sites

Diatoms

Only 2 sites have diatoms so not doing anything with this at the moment.

## # A tibble: 2 × 2
##   site            n
##   <fct>       <int>
## 1 Ruckle          6
## 2 Hailstorm 2     2

Functional groupings

Grazed Algae

exp_grazed_poisson <- glmmTMB(round(grazed) ~ treatment * region +  (1|region/site), data = exp_data, family = poisson(link = "log"))

exp_grazed_NB <- update(exp_grazed_poisson, family = nbinom2) 

exp_grazed_NB1 <- update(exp_grazed_poisson, family = nbinom1) 

exp_grazed_gaussian <- glmmTMB((grazed+1) ~ treatment * region +  (1|region/site), data = exp_data, family = gaussian(link = "log"))

exp_grazed_tweedie <- glmmTMB(grazed ~ treatment * region +  (1|region/site), data = exp_data, family = tweedie(link = "log")) 

# models with zero-inflation
exp_grazed_zipoisson <- glmmTMB(round(grazed) ~ treatment * region +  (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))

exp_grazed_ziNB <- update(exp_grazed_zipoisson, family = nbinom2) 

exp_grazed_ziNB1 <- update(exp_grazed_zipoisson, family = nbinom1) # best model

Checking assumptions:

##                      dAIC   df
## exp_grazed_ziNB1        0.0 8 
## exp_grazed_tweedie     11.0 8 
## exp_grazed_NB1         18.7 7 
## exp_grazed_ziNB        41.7 8 
## exp_grazed_NB          81.0 7 
## exp_grazed_gaussian   135.6 7 
## exp_grazed_zipoisson  542.0 7 
## exp_grazed_poisson   1291.1 6

Adding a dispersion formula for treatment*site

Model outputs:
coefficient estimates for model fitted to Grazed Algae field experiment data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.0771477 0.1244534 32.7604248 0.0000000
treatmentexclusion -0.1053661 0.1454185 -0.7245718 0.4687148
regionHigh -1.2498128 1.0176401 -1.2281482 0.2193914
treatmentexclusion:regionHigh 1.6313215 1.0190669 1.6007992 0.1094214
p values for model fitted to Grazed Algae field experiment data
term statistic df p.value
(Intercept) 1073.2454346 1 0.0000000
treatment 0.5250043 1 0.4687148
region 1.5083479 1 0.2193914
treatment:region 2.5625582 1 0.1094214
tukey adjusted p values for pairwise comparisons from model fitted to Grazed Algae field experiment data
contrast region estimate SE df t.ratio p.value
control - exclusion Low 0.1053661 0.1454185 65 0.7245718 0.4713150
control - exclusion High -1.5259554 1.0092699 65 -1.5119399 0.1353961

Nongrazed Algae

exp_nongrazed_poisson <- glmmTMB(round(nongrazed) ~ treatment * region +  (1|region/site), data = exp_data, family = poisson(link = "log"))

exp_nongrazed_NB <- update(exp_nongrazed_poisson, family = nbinom2) 

exp_nongrazed_NB1 <- update(exp_nongrazed_poisson, family = nbinom1) 

exp_nongrazed_gaussian <- glmmTMB((nongrazed+1) ~ treatment * region +  (1|region/site), data = exp_data, family = gaussian(link = "log"))

exp_nongrazed_tweedie <- glmmTMB(nongrazed ~ treatment * region +  (1|region/site), data = exp_data, family = tweedie(link = "log")) # best model 

# models with zero-inflation
exp_nongrazed_zipoisson <- glmmTMB(round(nongrazed) ~ treatment * region +  (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))

exp_nongrazed_ziNB <- update(exp_nongrazed_zipoisson, family = nbinom2) 

exp_nongrazed_ziNB1 <- update(exp_nongrazed_zipoisson, family = nbinom1)

Checking model assumptions:

##                         dAIC  df
## exp_nongrazed_tweedie     0.0 8 
## exp_nongrazed_ziNB1       6.8 8 
## exp_nongrazed_NB1         7.2 7 
## exp_nongrazed_ziNB       13.4 8 
## exp_nongrazed_NB         14.9 7 
## exp_nongrazed_zipoisson  15.7 7 
## exp_nongrazed_poisson    54.4 6 
## exp_nongrazed_gaussian  142.6 7

coefficient estimates for model fitted to Non-grazed Algae field experiment data
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.4324789 0.5585112 0.7743425 0.4387283
treatmentexclusion -0.4258574 0.3534379 -1.2049000 0.2282419
regionHigh -0.7537305 0.8135738 -0.9264440 0.3542153
treatmentexclusion:regionHigh -0.4247930 0.6973658 -0.6091394 0.5424320
p values for model fitted to Non-grazed Algae field experiment data
term statistic df p.value
(Intercept) 0.5996063 1 0.4387283
treatment 1.4517840 1 0.2282419
region 0.8582985 1 0.3542153
treatment:region 0.3710508 1 0.5424320
tukey adjusted p values for pairwise comparisons from model fitted to Non-grazed Algae field experiment data
contrast region estimate SE df t.ratio p.value
control - exclusion Low 0.4258574 0.3534379 76 1.204900 0.2319805
control - exclusion High 0.8506504 0.6012904 76 1.414708 0.1612383

Summary for all models:

Model descriptions for survey data
species distribution_family zero_inflated dispersion_formula month region month_region
limpets negative binomial 1 yes n/a yes yes yes
littorines negative binomial 1 no salinity*month no no no
balanus negative binomial 1 no salinity*month yes yes yes
chthamalus negative binomial 2 no salinity yes yes yes
red algae negative binomial 1 no site no yes no
green algae negative binomial 1 no n/a yes yes yes
brown algae negative binomial 1 no site yes no no
grazed negative binomial 2 no n/a yes no yes
nongrazed negative binomial 1 no site yes no no
Model descriptions for survey data
species distribution_family zero_inflated dispersion_formula treatment region treatment_region
balanus tweedie no treatment*site yes yes (site) yes
chthamalus negative binomial 1 yes n/a yes yes yes
all algae negative binomial 1 yes region yes yes yes
red algae tweedie no n/a no no no
green algae negative binomial 1 yes region no yes no
brown algae tweedie no n/a no no no
grazed negative binomial 1 yes treatment*site no yes no
nongrazed tweedie no n/a no no no